##### Data Fig5 - Geographic variation of mutagenic exposures in kidney cancer genomes
library(RColorBrewer)
library(pheatmap)
library(tidyverse)
library(viridis)
library(reshape2)
filter = dplyr::filter

#Load List of Drivers
drivers <- read.csv(file="Fig5_Driver_Input.csv")

#Separate Drivers and create short version
drivers_summary <- as.data.frame(cbind(drivers$Sample, drivers$Gene))
colnames(drivers_summary) = c("SampleID", "Gene")
#identify duplicated
duplciates <- as.data.frame(drivers_summary[duplicated(drivers_summary[,1:2]),])
#remove duplicates
drivers_summary_unique <- drivers_summary %>%
  distinct(SampleID, Gene, .keep_all = TRUE)

#####Fig5a
#calculate overall frequency
driver_frequency <- as.data.frame(table(drivers_summary_unique$Gene))
colnames(driver_frequency)[1] <- "Gene"
driver_frequency$Percent = driver_frequency$Freq / 962 *100

#Restrict to most frequent genes
driver_frequency_subset <- filter(driver_frequency, Freq > 10)

#Plot
pdf(file="Fig5a.pdf", width = 4, height = 2.5)
ggplot(driver_frequency_subset, aes(y=Percent, x=reorder(Gene, -Percent))) +
  geom_col(aes(fill = reorder(Gene, -Percent)), colour="black", position = "dodge") +
  scale_fill_viridis(discrete = TRUE) +
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  scale_y_continuous(limits = c(0,80), expand = c(0, 0)) +
  labs(y="% of cases") + 
  theme(axis.title.x = element_blank()) +
  theme(axis.title = element_text(size = 8)) +
  theme(axis.text = element_text(size = 8)) +
  theme(axis.text.x = element_text(angle=90)) +
  theme(legend.position = "none") 
dev.off()

#####Fig5b
#Load Country Data
Country <- read.csv(file="Fig5_Country_Input.csv")
Country <- as.data.frame(lapply(Country, gsub, pattern = "United Kingdom", replacement = "UK", fixed = TRUE))
Country <- as.data.frame(lapply(Country, gsub, pattern = "Czech Republic", replacement = "Czechia", fixed = TRUE))

#Removing countries with less than 20 cases
heatmap_data = Country
heatmap_data <- filter(heatmap_data, country != "Thailand")
heatmap_data <- filter(heatmap_data, country != "Poland")
heatmap_data <- filter(heatmap_data, country != "Lithuania")

#Make Heatmap Data
VHL_cases <- filter(drivers_summary_unique,Gene == "VHL")
SETD2_cases <- filter(drivers_summary_unique,Gene == "SETD2")
BAP1_cases <- filter(drivers_summary_unique,Gene == "BAP1")
PBRM1_cases <- filter(drivers_summary_unique,Gene == "PBRM1")

heatmap_data$VHL <- ifelse(heatmap_data$donor_id %in% VHL_cases$SampleID, 1, 0)
heatmap_data$PBRM1 <- ifelse(heatmap_data$donor_id %in% PBRM1_cases$SampleID, 1, 0)
heatmap_data$SETD2 <- ifelse(heatmap_data$donor_id %in% SETD2_cases$SampleID, 1, 0)
heatmap_data$BAP1 <- ifelse(heatmap_data$donor_id %in% BAP1_cases$SampleID, 1, 0)

#Summarise and plot heatmap
summary <- heatmap_data %>% group_by(country) %>% summarise_if(is.numeric, mean)
summary <- tibble::column_to_rownames(summary, var="country")
summary <- summary*100

pdf("Fig5b.pdf")
pheatmap(summary, 
         cluster_rows = FALSE, 
         cluster_cols = FALSE, 
         border_color = "black", 
         cellwidth = 15,
         cellheight = 15,
         angle_col = 45,
         main = "% of cases",
         display_numbers = round(summary,0))
dev.off()

#####Fig5c
input <- read.csv(file="Fig5_Spectra_Input.csv")

#Spectra Plot Function
plot_driver_spectra = function(input, plot_title) {
  #Collect Data
  mutations = input
  #Define contexts
  matrix_contexts <- c('A[C>A]A', 'A[C>A]C', 'A[C>A]G', 'A[C>A]T', 'C[C>A]A', 'C[C>A]C', 'C[C>A]G', 'C[C>A]T',
                       'G[C>A]A', 'G[C>A]C', 'G[C>A]G', 'G[C>A]T', 'T[C>A]A', 'T[C>A]C', 'T[C>A]G', 'T[C>A]T',
                       'A[C>G]A', 'A[C>G]C', 'A[C>G]G', 'A[C>G]T', 'C[C>G]A', 'C[C>G]C', 'C[C>G]G', 'C[C>G]T', 
                       'G[C>G]A', 'G[C>G]C', 'G[C>G]G', 'G[C>G]T', 'T[C>G]A', 'T[C>G]C', 'T[C>G]G', 'T[C>G]T',
                       'A[C>T]A', 'A[C>T]C', 'A[C>T]G', 'A[C>T]T', 'C[C>T]A', 'C[C>T]C', 'C[C>T]G', 'C[C>T]T',
                       'G[C>T]A', 'G[C>T]C', 'G[C>T]G', 'G[C>T]T', 'T[C>T]A', 'T[C>T]C', 'T[C>T]G', 'T[C>T]T',
                       'A[T>A]A', 'A[T>A]C', 'A[T>A]G', 'A[T>A]T', 'C[T>A]A', 'C[T>A]C', 'C[T>A]G', 'C[T>A]T',
                       'G[T>A]A', 'G[T>A]C', 'G[T>A]G', 'G[T>A]T', 'T[T>A]A', 'T[T>A]C', 'T[T>A]G', 'T[T>A]T',
                       'A[T>C]A', 'A[T>C]C', 'A[T>C]G', 'A[T>C]T', 'C[T>C]A', 'C[T>C]C', 'C[T>C]G', 'C[T>C]T',
                       'G[T>C]A', 'G[T>C]C', 'G[T>C]G', 'G[T>C]T', 'T[T>C]A', 'T[T>C]C', 'T[T>C]G', 'T[T>C]T',
                       'A[T>G]A', 'A[T>G]C', 'A[T>G]G', 'A[T>G]T', 'C[T>G]A', 'C[T>G]C', 'C[T>G]G', 'C[T>G]T',
                       'G[T>G]A', 'G[T>G]C', 'G[T>G]G', 'G[T>G]T', 'T[T>G]A', 'T[T>G]C', 'T[T>G]G', 'T[T>G]T')
  #Count Occurances
  freq=c(table(mutations$context))
  mutations_ctx=data.frame(ctxt=matrix_contexts,counts=freq[matrix_contexts])
  mutations_ctx <- tibble::column_to_rownames(mutations_ctx,"ctxt")
  mutations_ctx[is.na(mutations_ctx)] = 0
  rownames(mutations_ctx) <- NULL
  
  #Plot!
  # Specify Context Type
  sig_cat = c("C>A","C>G","C>T","T>A","T>C","T>G")
  ctx_vec = paste(rep(c("A","C","G","T"),each=4),rep(c("A","C","G","T"),times=4),sep="-")
  full_vec = paste(rep(sig_cat,each=16),rep(ctx_vec,times=6),sep=",")
  snv_context = paste(substr(full_vec,5,5), substr(full_vec,1,1), substr(full_vec,7,7), sep="")
  # Specify Vectors for plot colours
  col_vec_num <- rep(16,6)
  col_vec = rep(c("dodgerblue","black","red","grey70","olivedrab3","plum2"),each=16)
  # Convert to matrix
  sig_plot <- as.matrix(mutations_ctx)
  # Set up Signature Names and title
  sig_title <- colnames(sig_plot)
  # Set up Counts
  sample_counts <- colSums(mutations_ctx)
  # Set par settings
  par(xaxs='i', cex = 1, xpd = FALSE)
  #  define maxy
  max_prob <- sig_plot[,1];maxy = max(max_prob)
  #  call empty plot
  b = barplot(sig_plot[,1], col = NA, border = NA, axes = FALSE, las = 2, ylim=c(0,1.5*maxy), 
              names.arg = snv_context, cex.lab = 1.3, cex.names = 0.70, cex.axis = 2, 
              space = 1, ylab = "Mutation Count")  
  # add axis
  axis(2, at = pretty(0:(1.5*maxy), n = 3), col = 'grey90', las = 1, cex.axis = 1.5)
  # call columns
  b = barplot(sig_plot[,1], axes = FALSE, col=col_vec, add = T, border = NA, space = 1)
  # add box surronding plot
  box(lty = 1, col = 'grey90')
  # add title
  title(plot_title, line = -1.5, adj = 0.01, cex.main = 2)
  title(paste0("Total SNV: ", sample_counts), line = -1, adj = 0.99, cex.main = 1)
  # add rectangles and annotations on top of the plot
  par(xpd = NA)
  for (j in 1:length(sig_cat)) {
    xpos = b[c(sum(col_vec_num[1:j])-col_vec_num[j]+1,sum(col_vec_num[1:j]))]
    rect(xpos[1]-0.5, maxy*1.5, xpos[2]+0.5, maxy*1.6, border=NA, col=unique(col_vec)[j])
    text(x=mean(xpos), pos=3, y=maxy*1.6, label=sig_cat[j], cex = 1.25)
  } 
}

#Subset only for cases with >10% COSMIC attribution to either AA signatures
AA_exposed <- read.csv(file="Fig5_AA_Input.csv")
input_AA <- subset(input, (sampleID %in% AA_exposed$Sample))
input_nonAA <-subset(input, (!sampleID %in% AA_exposed$Sample))
#Plot driver spectra in AA vs nonAA
pdf(file="Fig5c_AA.pdf", width=12,height=4)
plot_driver_spectra(input_AA,"AA exposed")
dev.off()
pdf(file="Fig5c_Unexposed.pdf", width=12,height=4)
plot_driver_spectra(input_nonAA,"Unexposed")
dev.off()


#####Fig5d
#Plot Difference in T>A as a proportion
input_TtoA <- filter(input, sub == "T>A")
input_TtoA$group <- ifelse(input_TtoA$sampleID %in% AA_exposed$Sample, "Exposed", "Unexposed")
input_TtoA_table <- as.data.frame(table(input_TtoA$group))
input_TtoA_table$no_variants = c(nrow(input_AA), nrow(input_nonAA))
input_TtoA_table$percentage <- input_TtoA_table$Freq / input_TtoA_table$no_variants * 100
input_TtoA_table <- cbind(input_TtoA_table[,c(1,4)])
colnames(input_TtoA_table) = c("Status", "T>A Mutations")
#Prepare to Plot
input_TtoA_melt <- melt(input_TtoA_table)
input_TtoA_melt$Status <- factor(input_TtoA_melt$Status, levels = c("Exposed","Unexposed"))
#Plot
pdf(file="Fig5d.pdf", height = 4, width =4)
ggplot(input_TtoA_melt,aes(y=value,x=variable)) +
  geom_col(aes(fill = Status), colour="black", position = "dodge") +
  scale_fill_viridis(discrete = TRUE) +
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  labs(y="% of total mutations") + 
  theme(axis.title.x = element_blank()) +
  theme(axis.title = element_text(size = 12)) +
  theme(axis.text = element_text(size = 12)) +
  theme(legend.position = "bottom") +
  theme(legend.title = element_blank()) +
  theme(legend.text=element_text(size=12)) 
dev.off()

#####Fig5e
#filter for VHL drivers only
input_VHL <- filter(input, gene == "VHL")
#Subset only for cases with >10% COSMIC attribution to either AA signatures
input_VHL_AA <- subset(input_VHL, (sampleID %in% AA_exposed$Sample))
input_VHL_nonAA <-subset(input_VHL, (!sampleID %in% AA_exposed$Sample))
#Plot driver spectra in AA vs nonAA
pdf(file="Fig5e_AA.pdf", width=12,height=4)
plot_driver_spectra(input_VHL_AA,"AA exposed")
dev.off()
pdf(file="Fig5e_Unexposed.pdf", width=12,height=4)
plot_driver_spectra(input_VHL_nonAA,"Unexposed")
dev.off()

#####Fig5f
#Plot Difference in T>A as a proportion
input_TtoA_VHL <- filter(input_VHL, sub == "T>A")
input_TtoA_VHL$group <- ifelse(input_TtoA_VHL$sampleID %in% AA_exposed$Sample, "Exposed", "Unexposed")
input_TtoA_VHL_table <- as.data.frame(table(input_TtoA_VHL$group))
input_TtoA_VHL_table$no_variants = c(nrow(input_VHL_AA), nrow(input_VHL_nonAA))
input_TtoA_VHL_table$percentage <- input_TtoA_VHL_table$Freq / input_TtoA_VHL_table$no_variants * 100
input_TtoA_VHL_table <- cbind(input_TtoA_VHL_table[,c(1,4)])
colnames(input_TtoA_VHL_table) = c("Status", "T>A Mutations")
#Prepare to Plot
input_TtoA_VHL_melt <- melt(input_TtoA_VHL_table)
input_TtoA_VHL_melt$Status <- factor(input_TtoA_VHL_melt$Status, levels = c("Exposed","Unexposed"))
#Plot
pdf(file="Fig5f.pdf", height = 4, width =4)
ggplot(input_TtoA_VHL_melt,aes(y=value,x=variable)) +
  geom_col(aes(fill = Status), colour="black", position = "dodge") +
  scale_fill_viridis(discrete = TRUE) +
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  labs(y="% of total mutations") + 
  theme(axis.title.x = element_blank()) +
  theme(axis.title = element_text(size = 12)) +
  theme(axis.text = element_text(size = 12)) +
  theme(legend.position = "bottom") +
  theme(legend.title = element_blank()) +
  theme(legend.text=element_text(size=12)) 
dev.off()







